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Abstract. We study the chemostat model for one species competing 
for one nutrient using a Lyapunov-type analysis. We design the dilution 
rate function so that all solutions of the chemostat converge to a pre- 
scribed periodic solution. In terms of chemostat biology, this means that 
no matter what positive initial levels for the species concentration and 
nutrient are selected, the long term species concentration and substrate 
levels closely approximate a prescribed oscillatory behavior. This is sig- 
nificant because it reproduces the realistic ecological situation where the 
species and substrate concentrations oscillate. We show that the stabil- 
ity is maintained when the model is augmented by additional species 
that are being driven to extinction. We also give an input-to-state sta- 
bility result for the chemostat tracking equations for cases where there 
are small perturbations acting on the dilution rate and initial concentra- 
tion. This means that the long term species concentration and substrate 
behavior enjoys a highly desirable robustness property, since it contin- 
ues to approximate the prescribed oscillation up to a small error when 
there are small unexpected changes in the dilution rate function. 
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1. Introduction. The chemostat model provides the foundation for much 
of current research in bio-engineering, ecology, and population biology [6, 
8, 10, 11, 21]. In the engineering literature, the chemostat is known as 
the continuously stirred tank reactor. It has been used for modeling the 
dynamics of interacting organisms in waste water treatment plants, lakes and 
oceans. In its basic setting, it describes the dynamics of species competing 
for one or more limiting nutrients. If there are n species with concentrations 
Xi for i = 1, . . . , n and just one limiting nutrient with concentration S and 
dilution rate D > 0, then the model takes the form 

S = D{S in - S) -^2(j, i (S)x i /'Yi 
±i = Xi(fj,i(S) - D), i = l,...,n 

where fii denotes the per capita growth rate of species i and p is the time 
derivative of any variable p. (In much of the paper, we simplify our notation 
by omitting the arguments of the functions. For instance, when no confusion 
can arise from the context, we denote S(t) simply by S.) The functions 
depend only on the nutrient concentration, and are zero at zero, continuously 
differentiable and strictly increasing, although non-monotone functions have 
been the subject of research as well. The conversion of nutrient into new 
biomass for each species i happens with a certain yield ji £ (0, 1) and the 
natural control variables in this model are the input nutrient concentration 
Si n and the dilution rate D. The latter variable is defined as the ratio of 
the volumetric flow rate F (with units of volume over time) and the reactor 
volume V r which is kept constant. Therefore it is proportional to the speed 
of the pump that supplies the reactor with fresh medium containing the 
nutrient. The equations (1) are then straightforwardly obtained from writing 
the mass-balance equations for the total amounts of the nutrient and each 
of the species, assuming the reactor content is well-mixed. The full model 
(1) is illustrated in Figure 1. 
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Figure 1. Chemostat 
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In the present work, we consider the case where there is just one species 
with concentration x, in which case the equations (1) take the form 

f S = D(S in - S) - ii(S)x/rf , 2 s 
\ x = x(n{S) -D) 

(but see Theorem 3 below for results on chemostats with disturbances, and 
Section 8 for models involving several species). We assume Si n is a given 
positive constant, while the per capita growth rate \i is a Monod function 
(which is also known as a Michaelis-Menten function) taking the form 

m = (3) 

a + b 

for certain positive constants m and a that we specify later. The dilution 
rate is an appropriate continuous positive periodic function we also specify 
below. Since S > when S > is near zero, one can readily check that (2) 
leaves the domain of interest 

X := (0,oo) x (0,oo) 

positively invariant (i.e., trajectories for (2) starting in X remain in X for 
all future times); see Theorem 3 for a more general invariance result for 
perturbed chemostats. 

Since we are taking Si n to be a fixed positive constant, we rescale the 
variables to reduce the number of parameters. Using the change of variables 

S x 

S =-q-l -q ) U(S)=n(S in S) 

and dropping bars, we eliminate Si n and 7 and so obtain the new dynamics 

f S = D(l-S)- »(S)x 
\ x = x(n(S)-D) 

again evolving on the state space X = (0, oo) 2 . Motivated by the realistic 
ecological situation where the species concentrations are oscillating, we solve 
the following biological problem: 

Biological Problem Bl: For a prescribed oscillatory behavior for 
the species concentration and substrate level given by a time-periodic 
pair (S r (t),x r (t)), design a dilution rate function D(t) such that if 
this choice of D(t) is used in the chemostat (4), then all solution pairs 
(S(t),x(t)) for the substrate levels and corresponding species levels 
obtained from solving (4) (i.e. for all possible initial values) closely 
approximate (S r (t), x r (t)) for large times t. 

See also Problem (SP) in Section 4 below for a precise mathematical state- 
ment of the preceding problem. In the language of control theory, solving 
Biological Problem Bl means we will prove the stability of a suitable peri- 
odic reference signal for the species concentration t \— > x r (t) in (4) using an 
appropriate time-periodic dilution rate D(t); see [15] for the fundamental 
ideas from control theory we need in the sequel. 
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Since D{t) is proportional to the speed of the pump which supplies the 
chemostat with medium containing nutrient, implementation of the pre- 
scribed oscillatory behavior requires that we control the pump in a very 
precise way. In practice this control process is prone to errors, and the ac- 
tual pump speed will be subject to small fluctuations which we will model 
by replacing D(t) by D (t) + ui(t) in the chemostat equations, where u\(t) 
models the error. It is therefore of interest to study the effect of these small 
fluctuations on the periodic behavior. Preferably this effect will be small, 
and the resulting behavior is not too different from the prescribed periodic 
behavior. We will show that this is indeed the case by actually quantifying 
how small the deviations are, relying on the well-known control-theoretic no- 
tion of Input-to-State Stability or ISS; see [23, 24] and Remark 2 for details 
about ISS. Summarizing this a bit more formally, we will solve the following 
biological problem (which we state in a more precise mathematical way in 
Section 5): 

Biological Problem B2: For the prescribed oscillatory behavior 
(S r (t),x r (t)) and dilution rate D(t) obtained in Biological Problem 
Bl, quantify how the substrate and species levels (S(t),x(t)) in the 
chemostat model (4) are affected by unexpected changes in the dilution 
rate, and show that the convergence of (S(t), x(t)) to the oscillatory 
behavior (S r (t),x r (t)) is robust to small changes in D(t). 

Our solution to Biological Problem B2 will be a special case of our more 
general input-to-state stability result for the chemostat tracking equations, 
assuming the dilution rate and initial concentration are both perturbed by 
noise terms of small enough magnitude. 

In the next section, we briefly review the literature focusing on what 
makes our approach different. In Section 3, we fix the reference signal we 
wish to track. In Section 4, we precisely formulate the definitions and the 
stability problem we are solving. We state our main stability theorem in 
Section 5 and we discuss the significance of our theorem in Section 6. We 
prove our stability result in Section 7. In Section 8, we show that the stability 
is maintained when there are additional species that are being driven to 
extinction. We validate our results in Section 9 using a numerical example. 
We conclude in Section 10 by summarizing our findings. 

2. Review of the Literature and Comparison with Our Results. 

The behavior of the system (1) is well understood when Si n and D are 
positive constants, as well as cases where n = 2 and either of these control 
variables is held fixed while the other is periodically time-varying. See [13, 
20] for periodic variation of Si n and [4] for periodic variation of D and the 
general reference [21] on chemostats. When both Si n and D are constants, 
the so-called "competitive exclusion principle" holds, meaning that at most 
one species survives. Mathematically this translates into the statement that 
system (1) has a steady state with at most one nonzero species concentration, 
which attracts almost all solutions; see [21]. This result has triggered much 
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research to explain the discrepancy between the (theoretical) competitive 
exclusion principle and the observation that in real ecological systems, many 
species coexist. 

The results on the periodically-varying chemostat mentioned above should 
be seen as attempts to explain this paradox. They involve chemostats with 
n = 2 species, and their purpose is to show that an appropriate periodic 
forcing for either Si n (t) or D(t) can make the species coexist, usually in the 
form of a (positive) periodic solution. Few results on coexistence of n > 2 
species are available. An exception is [19], where a periodic function Si n (t) 
is designed (with D kept fixed) so that the resulting system has a (positive) 
periodic solution with an arbitrary number of coexisting periodically varying 
species. The stability properties of this solution are not known. 

More recent work has explored the use of state-dependent but time invari- 
ant feedback control of the dilution rate D to generate coexistence; see [6, 8] 
for monotone growth rate functions in the n = 2 species case, and [7] for 
the n = 3 species case. The paper [11] considers feedback control when the 
growth rate functions are non-monotone. In [12], [16], and [18], coexistence 
is proved for models taking into account intra-specific competition. In these 
models, the usual growth functions fJ.i(S) are replaced by functions Hi(S, Xi) 
which are decreasing with respect to the variable X{. All the results discussed 
so far apply to a more general model than (4) involving n > 1 species. This 
is because the main purpose of these papers is to investigate environmental 
conditions under which the competitive exclusion principle fails and several 
species can coexist. 

Here we will not consider any coexistence problems. Our main objective 
is to provide a proof of stability of a periodic solution based on a Lyapunov- 
type analysis and to investigate the robustness properties of the periodic 
solution with respect to perturbations. As an illustration we show that the 
stability of the periodic solution is robust with respect to additional species 
that are being driven to extinction, or to small disturbances on the initial 
nutrient concentration or dilution rate. These features set our work apart 
from the known results on periodically forced chemostat models which do 
not rely on the construction of a Lyapunov function. Proving stability in 
the chemostat usually relies on reduction and monotonicity arguments, and 
not so often on Lyapunov functions (but see for instance Theorem 4.1 in [21] 
which uses a Lyapunov function introduced in [14] and more recently [12]). 

Finally we point out that closely related to our results is [10] where a 
single-species chemostat with a continuous and bounded (but otherwise ar- 
bitrary) function Si n (t) and constant dilution rate is investigated; there it is 
shown that two positive solutions converge to each other. However, the proof 
is not based on a Lyapunov function. The advantage of having a Lyapunov 
function is that it can be used to quantify the effect of additional noise terms 
on the stability of the unperturbed dynamics. In fact, to our knowledge, our 
work provides the first input-to-state stability analysis of chemostats whose 
dilution rates and initial concentrations are perturbed by small noise; see 
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Remark 2 for a discussion on the importance of input-to-state stability in 
control theory and engineering applications. 



3. Choosing a Reference Trajectory. We first choose the dilution rate 
D = D(t) that will give a reference trajectory (S r (t),x r (t)) for (4) which 
we show to be stable. We assume a growth rate with constants m,a > 
as follows, in which the constants a and m and the variable S are all 
dimensionless by the change of coordinates used to obtain the normalized 
equations (4) so the units do not matter: 



mS 
a + S' 



where m > 4a + 1 . 



(5) 



For the sake of computational simplicity, we choose a sinusoidal reference 
trajectory but the extension to more general reference trajectories can be 
handled by similar methods; see Remark 5 for details. Simple calculations 
show that (4) admits the trajectory 



(S r {t),x r (t)) = (\-\ cos (*)> \ + \ cos 0) 
which we refer to as a reference trajectory when we choose 



D(t) 



x r {t) 



+ H(l-x T {t)) 



sin(t) 



m (2 — cos(f)) 



x r (t) 2 + cos(t) 4a + 2-cos(£)' 

Condition (5) then provides constants D,D a > such that 

D < D{t) < D 

for all t > 0: 



(6) 



(7) 



D = l + 



3m 



and D r , 



m 



1 . 



4a + 3 4a + 1 

See Figure 2 for the graph of D(t) for m = 10 and a = \. 



(8) 




Figure 2. Dilution Rate D(t) for the Chemostat from (7) 



STABILITY IN THE PERTURBED CHEMOSTAT 



7 



4. Definitions and Statement of Stability Problem. We wish to solve 
the following stability problem which is merely a restatement of Biological 
Problem Bl above in precise control theoretic terms: 

(SP) Given any trajectory (S, x) : [0, oo) — ► (0, oo) 2 for (4) correspond- 
ing to the dilution rate D(t) from (7) and /i as in (5) (i.e. for 
any initial value for (S, x)), show that the corresponding deviation 
(S(t),x(t)) := (S(t) — S r (t),x(t)—x r (t)) of(S,x) from the reference 
trajectory (6) asymptotically approaches (0,0) as t — > +oo. 

We will solve (SP) by proving a far more general tracking result for a 
single species chemostat acted on by a disturbance vector u = (u\,U2) ■ 
[0, oo) — ► R 2 as follows: 

S(t) = [D(t) + Ul (t)](l + u 2 (t)-S(t))-»(S(t))x(t) 
x(t) = x(t)[fi(S(t))-D(t)- Ul (t)} 

We will quantify the extent to which the reference trajectory (6) tracks the 
trajectories of (9) which will solve Biological Problem B2 from the intro- 
duction. To this end, we need to introduce a priori bounds on u\ and ui\ 
see Remark 8. Our main theoretical tool will be the input-to-state stability 
(ISS) property [22] which is one of the central paradigms of current research 
in nonlinear stability analysis; see Remark 2. The relevant definitions are 
as follows. 

We let /Coo denote the set of all continuous functions 7 : [0, 00) — > [0, 00) 
for which (i) 7(0) = and (ii) 7 is strictly increasing and unbounded. We let 
ICC denote the class of all continuous functions (3 : [0, 00) x [0, 00) — ► [0, 00) 
for which (I) (3(-,t) G /Coo for each t > 0, (II) 0(s, •) is non-increasing for 
each s > 0, and (III) (3(s, t) — > as t — > +00 for each s > 0. Consider a 
general control-affine dynamic 

y = F(y,t) + G(y,t)u, y£0, uGU (10) 

evolving on a given subset O C W 1 where U is a given subset of Euclidean 
space. (Later we specialize to dynamics for the chemostat.) For each t a > 
and y G O, let y(t; t Q , y , a) denote the solution of (10) satisfying y(t ) = y 
for a given control function a G U := {measurable essentially bounded a : 
[0,oo) — > U}; i.e. the solution of the initial value problem 

y(t) = F{y(t),t) + G(y{t),t)a(t) a.e. t>t a , y{t ) = y Q . 

We always assume that such solutions are uniquely defined on all of [t , 00) 
(i.e., (10) is forward complete and O is positively invariant for this sys- 
tem) and that there exists O G /Coo such that \F(y,t)\ + \G(y, t)\ < &(\y\) 
everywhere, where | • | is the usual Euclidean norm. For example, 

[t , 00) 3 1 1 ^ (S(t; t , (S , x ), a),x(t; t Q , (S Q , x Q ), a)) 

is the solution of (9) for the disturbance u = (u\,U2) = a(t) satisfying the 
initial condition (S(t ),x(t )) = (S ,x ). 
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Definition 1. We call (10) input-to-state stable (ISS) provided there exist 
(3 G ICC and 7 G /Coo such that 



for all t > t Q , t Q > 0, Do G 0, and a €U. 

Here |a|oo denotes the essential supremum of a G W. By causality, the ISS 
condition (11) is unchanged if |a|oo is replaced by the essential supremum 
l a l[t ,t] °f a restricted to [t Q ,t]. In particular, (11) says y(t;t ,y , Z) — ► 
as t — ► +00 for all initial values y Q and initial times t , where Z is the zero 
disturbance a{t) = 0. 

Remark 2. T/ie theory of ISS systems originated in [22]. 755 theory pro- 
vides the foundation for much current research in robustness analysis and 
controller design for nonlinear systems, and has also been used extensively in 
engineering and other applications [1, 2, 5, 17, 22, 24]. The ISS approach can 
be viewed as a unification of the operator approach of Zames (e.g. [25, 26] J 
and the Lyapunov state space approach. The operator approach involves 
studying the mapping (t ,y ,a) 1— > y(-;t ,y ,a) of initial data and control 
functions into appropriate spaces of trajectories, and it has the advantages 
that it allows the use of Hilbert or Banach space techniques to generalize 
many properties of linear systems to nonlinear dynamics. By contrast, the 
state space approach is well suited to nonlinear dynamics and lends itself to 
the use of topological or geometric ideas. The ISS framework has the advan- 
tages of both of these approaches including an equivalent characterization 
in terms of the existence of suitable Lyapunov-like functions; see Remark 7 
below. For a comprehensive survey on many recent advances in ISS theory 
including its extension to systems with outputs, see [24]. 

To specify the bound u on our disturbances u = (111,112), we use the 
following constants whose formulas will be justified by the proof of our main 
stability result: 



5. Statement of Theorem. From now on, we assume the disturbance 
vector u = (7/1,7/2) in (9) takes all of its values in a fixed square control set 
of the form 



where c, k, and C\ are in (12) (but see Remark 8 for related results under less 
stringent conditions on the disturbance values). We will prove the following 
robustness result: 



y(t;t ,y ,a)\ < f3(\y \,t - t Q ) + 7(Moo) 



(11) 




(12) 



U := {(7/1,7/2) G M 2 : \ui\ < u, 1 7/2 1 < u} where 




(13) 
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Theorem 3. Choose D{t), fi, and (S r ,x r ) as in (5)-(7). Then the corre- 
sponding solutions of (9) satisfy 

\S o >0 & x o >0 & t>t o >0 & a€U] 

~ J (14) 

[ S{t;t , (S ,x ),a) > & x{t;t , (S ,x ),a) > ] . 

Moreover, there exist (3 G XX and 7 G /Coo sitc/i i/iai £/ie corresponding 
transformed error vector 

y{t;t ,y ,a) := 

{S{t; t , (So, so), «) - SV(i), ln(x(t; i , (So, x ),a)) - \n(x r {t))) 

satisfies the ISS estimate (11) for all a £ U , to > 0, t > to, So > 0, and 
x > 0, where y a = {S ,x ). 



6. Discussion on Theorem 3. Before proving the theorem, we discuss 
the motivations for its assumptions, and we interpret its conclusions from 
both the control theoretic and biological viewpoints. 

Remark 4. Condition (14) says (0,oo) 2 is positively invariant for (9). One 
may also prove that [0, oo) 2 is positively invariant for (9), as follows. Sup- 
pose the contrary. Fix t Q > 0, x Q > 0, S > 0, and a G U for which the cor- 
responding trajectory (S(t),x(t)) for (9) satisfying {S{t ),x{t )) = (S ,x ) 
exits [0, oo) 2 in finite time. This provides a finite constant t\ := max{£ > 
t a : (S(t),x(t)) G [0,oo) 2 Vi G [t , t]}-. Then S{t±) = 0, since otherwise 
S(t±) > and x{t) = for all t > t\ and then we could use the continuity of 
S to contradict the maximality oft\. Since u < min{l, D Q /2}, the continuity 
of S and x and the fact that S{t±) = provide a constant e > such that 
1 + u 2 {t) - S(t) > (1 - u)/2 and n{S{t))x{t) < D (l - u)/8 for (almost) 
all t G [*i,ti + e], hence also S{t) > D (l - u)/8 > for all t G [ti,h + e] 
(since D{t) + u\{t) > D Q /2 everywhere). Hence, S{t) > S{t±) = for 
all t G [ti,t\ + e]. Since x{t) clearly stays in [0, 00), this contradicts the 
maximality oft±. The positive invariance of [0, oo) 2 follows. 

Remark 5. Theorem 3 says that in terms of the error signals y, any com- 
ponentwise positive trajectory of the unperturbed chemostat dynamics (9) 
converges to the nominal trajectory (6), uniformly with respect to initial 
conditions. This corresponds to putting a = in (11). It also provides the 
additional desirable robustness property that for an arbitrary TJ-valued con- 
trol function a G U, the trajectories of the perturbed chemostat dynamics 
(9) are "not far" from (6) for large values of time. In other words, they 
"almost" track (6) with a small overflow 7(|a|oo) from the ISS inequality 
(11). Similar results can be shown for general choices of x r and D{t). For 
example, we can choose any x r {t) that admits a constant i > such that 

3 

max{^, |£ r (t)|} < x r {t) < - 
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for all t > and S r = 1 — x r . In this case, we take the dilution rate 

x r (t) 

which is again uniformly bounded above and below by positive constants. The 
proof of this more general result is similar to the proof of Theorem 3 we give 
below except with different choices of the constants c and k. 

Remark 6. The robustness result 

| {S{t; t , (Sq,x ), a) - S r {t),ln(x(t; t , (S , x ),a)) - ln(x r (t))) | 

(15) 

< 0(\(So,x o )\,t- t ) +7(|a|[t 0l *]) 
of Theorem 3 differs from the classical ISS condition in the following ways: 

1. For biological reasons, negative values of the nutrient level S and the 
species level x do not make physical sense. Hence, only componentwise 
positive solutions are of interest. Therefore, (15) is not valid for all 
(Sq,xo) S M 2 but rather only for (So,xo) G (0, oo) x (0, oo). 

2. Our condition (15) provides an estimate on the transformed error com- 
ponent ln(x(i;toi {Sq,xq),o)) — \n.(x r (t)) instead of the more standard 
error x(t; to, (So,xo),a) — x r (t) . Our reasons for using the transformed 
form of the error are as follows. The function ln(x) goes to —oo 
when x goes to zero. This property is relevant from a biological point 
of view. Indeed, in the study of biological systems, it is important 
to know if the concentration of the species is above a strictly posi- 
tive constant when the time is sufficiently large or if the concentra- 
tion admits zero in its omega limit set. In the first case, the species 
is called persistent. The persistency property is frequently desirable, 
and it is essential to know whether it is satisfied. Hence, the function 
ln(x(t;to,(So,xo),a)) — ln(x r (t)) has the desirable properties that (a) 
it goes to +oo if x(t; to, (So, xo), a) does, (b) it is equal to zero when 
x(t;to, (So,xq),u) is equal at time t to the value of x r , and (c) it goes 
to — oo if x(t; to, [So, xo), u) goes to zero. Therefore, roughly speaking, 
if the species faces extinction, then it warns us. 

Remark 7. Our proof of Theorem 3 is based on a Lyapunov type analysis. 
Recall that a C 1 function V : W 1 x [0, oo) -> [0, oo) is called an ISS Lyapunov 
function (ISS-LF) for (10) provided there exist 71,72,73,74 G K-oo such that 

1. 7i(M) < V(y,t)<j 2 (\y\) and 

2. V t (y,t) + V y (y,t)[F(y,t) + G(y,t)u] < - l3 {\y\) + 74(H) 

hold for all y £ O, t > 0, and u G U. The function V we will construct in 
the proof of Theorem 3 is not an ISS-LF for the chemostat error dynamics 
because of the specificities of the state space, which preclude the existence of 
the necessary functions 71,72 G /Coo in Condition 1 above. Hence, we can- 
not directly apply the result that the existence of an ISS Lyapunov function 
implies that the system is ISS e.g. [9, Theorem 1] to prove our theorem. 
Instead, we prove our Theorem 3 directly from the decay inequality satisfied 
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by the time derivative of V along the trajectories. The proof that the decay 
inequality implies ISS is very similar to that part of the proof of [9, Theorem 
1] , so we only sketch that part of our proof in the appendix. 

Remark 8. Our estimate (15) would not hold if we had instead chosen the 
full control set U = M. 2 . In fact, taking the disturbance a = (7/1,1*2) = 
(0,-1) and any initial condition (S(t ),x(t )) = (Sq,x ) G (0, oo) 2 in (9) 
would give S(t;to,(So,xo),a) — > and so also ln(x(t; to, (So, xo), <*))) — * 
—00 as t — ► +00 (since \ui(t)\ < u < D Q /2 < D(t)/2 almost everywhere). 
Therefore extinction would occur and (15) would not be satisfied. On the 
other hand, if our set U is replaced by Tj" := [—u, +u] 2 for any fixed constant 
u G (0, min{l, D Q }), then the chemostat error dynamics instead satisfies the 
less stringent integral ISS property; see Remark 9 below for details. 



7. Proof of Theorem 3. The proof of (14) is immediate from the structure 
of the dynamics (9) and the fact that u < 1 (which imply that S > when 
S > is sufficiently small); see Remark 4 for a similar argument. It remains 
to prove the ISS estimate (15) for suitable functions 3 G ICC and 7 G /Coo. 

Throughout the proof, all (in) equalities should be understood to hold 
globally unless otherwise indicated. Also, we repeatedly use the simple 
"(generalized) triangle inequality" relation 

pq < d P 2 + ^ 2 ( 16 ) 

for various choices of p > 0, q > 0, and d > that we specify later. 

Fix t a > 0, S > 0, x > 0, and a eU, and let [t Q , 00) B t (->• (S(t),x(t)) 
denote the corresponding solution of (9) satisfying (S(to),x(t )) = (S ,x ). 
For simplicity, we write a(t) as (7/1,1/2), omitting the time argument as 
before. We first write the error equation for the variables 

(2,£) = (z-Zr,S-£r) (17) 

where £ = lnx, z = S + x, z r (t) = S r (t) + x r (t) = 1, and £ r (t) = lnx r (t). 
One easily checks that 

z(t) = [D(t) + Ul (t)][l+u 2 (t)-z(t)] 
x(t) = x(t)[fi(z(t) - x(t)) - D(t) - ui(t)] . 

Therefore, since z r = 1 (which implies z r (t) = [D(t) + ui(t)][l — z r (t)]), our 
formula (5) for /z immediately implies that the (transformed) error 

(z,i)(t) = (z(t) - z r (t),lnx(t) - lnx r (t)) 

satisfies the chemostat error dynamics 

' z = ~[D(t) + Ul (t)][z - u 2 {t)} 

< , g-e^M^'-i) (18) 

J - ma (a + z _ e e)( a + Zr(t )_ e ^W) niW - 
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We are going to show that (18) has the Lyapunov function 

V(l z) = e L ^ - 1, where L 3 (f, z) = Li(f) + kL 2 (z), 
Li(0 = e« - £ - 1, and L 2 {z) = j^Z 2 

and k > is the constant defined in (12). From the explicit expressions 
z r (t) = 1 and e^ r ^ = \ + ^cos(t) > |, we deduce that the time derivative 
of L\ along trajectories of (18) satisfies 

Ll = ma r^ — £v - L — rrm-^- 1 " 1 ^) 

(a + z — e*)(a + z r (t) — e? r W) 

1 (J _1\2, 4 I £ _ I ||=| 

4a+2-cos(t) ^ c L > ~ 4a+2-cos(t) l c "MM , t 1N 

< ma — ft (et-l)ui(t) 

(a + 2 — J 

.(pi 1^2, 4mo I J _ i ||=| 



^ 4a+3^ / 1 4a+H " ' / f n \ /.\ 
^ (a + 2-e^) (e 1)UlW ' 

where we also used the fact that z — = S > 0. Since D(t) + u\(t) > D D — u 
everywhere, one readily checks that along the trajectories of (18), 

L 2 < -I 2 + c\l\\u 2 {t)\ < -\~z 2 + cu\(i) (20) 

where c = 2(D + u)(D — u)~ 1 , the constant c is defined by (12), and the last 
inequality used (16) with the choices p = \z\, q = \u 2 (t)\, and d = l/(2c). 
The fact that ^c 2 < c easily follows because u < \D . Along the trajectories 
of (18), 

i3 s ^ — !U_ (ee _ 1)Bi(t) (2i) 

— — KZ 2 + HCU^if). 

We distinguish between two cases. 

Case la : z(t) < 2. Then since z — = S > 0, we get 

• ma , ? so 4m . ; ... 

L 3 < -(e«-l 2 + e*-l z 

~ (4a + 3)(a + 2) V ; 4a + l' 11 1 ( 22 ) 

— (e^ — l)ni(t) — ^K2 2 + KCU 2 (t). 
Using the triangle inequality (16) with the choices 

D -ue_l| 0-151 an d d - ° (4fl + 1} 
p-|e 1|, and d - ^ + ^ + 2) , 

we deduce from (22) that 
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Case 2a : z(t) > 2. Since z(t)-e^ > 0, it follows that e~^z(t)-e^ > 0. 
Therefore, since x r > 1/4, 

e i(t) < e -^z{t) < e-^/^zit) =4z(t), m 
hence - 1 < - 1 < 4z{t) - I. 

Since z(t) > 2 and z r = 1, we have z(t) > 1. As z = z — z r = z — 1, condition 
(24) gives 

-z{t) < < 3 + 4z(t) < 7z(t), hence \e*® - 1| < 7z(t). (25) 

From this last inequality and the inequality z — > 0, we deduce from 
dropping the first term in (21) that 



• 28m n , c . K r> 



KCU 



< — L K ( e e_i)2 

200 v 7 



1 28m 



(26) 



2 



(e^ — l)ui + KCU2 ■ 



4 4a + 1 

We deduce from our choice (12) of n, (23), and (26) that in Cases la-2a, 

L 3 < -Ci[(e«-l) 2 + z 2 ]-(e«- l)m + kcu| (27) 
where Ci is defined in (12). 



Using (16) with p = |exp(£(i)) — 1|, q = \ui(t)\, and d = and then the 



upper bounds of u\ and tt2, we deduce from (27) that 

Ci 1 

L 3 < — — [(e $ - l) 2 + r 2 ] + — u|ui| + Kctt |u 2 1 

< -^"[(e^ — !) 2 + 5 2 ] -h where C 2 := + 2/ccJ it 

and where the last inequality used the relationship |u|i < 2|u| 2 between the 
1-norm and the 2-norm. We consider two additional cases. 
Case lb : (e*W - l) 2 + z 2 (t) > \. Then (28) gives 

V < eM^ (-^+C 2 \u(t)\\ . (29) 

Next notice that (13) and our choice of C\ G (0, 1] give 

u < hence V < -^-e L ^ < -Q-V&z). (30) 

Case 2b : (e^ - l) 2 + z 2 {t) < \. Then (f(t),5(i)) is in a suitable bounded 
set, so since 

F(L) = (e L - L - l)(e L - 1)~ 2 , G{L) = (e L - 1)L~ X 

are locally bounded when defined to be zero at L = 0, one can readily use 
(28) to compute constants C3 > and C4 > such that 

V < -C 3 L 3 (iz) + C 2 \u(t)\ < -C 4 V(£,z) + C 2 \u(t)\ (31) 
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where F was used to get C3 and G was used to get C4. It follows from 
(30)-(31) that, in Cases lb-2b, 

V < -C 5 V(i,z) + C 2 \u(t)\ (32) 

with C5 = inf {C4, g }• Condition (32) is the classical ISS Lyapunov func- 
tion decay condition for the transformed error dynamics evolving on our 
restricted state space. Therefore, a slight variant of the classical ISS argu- 
ments combined with (32) give the ISS estimate asserted by Theorem 3. For 
details, see the appendix below. This concludes the proof. 

Remark 9. If our control set U is replaced by the larger control set U" := 
[— u, +u] 2 for any fixed constant u £ (0, min{l, D }), then the error dynam- 
ics (18) instead satisfies the less stringent integral ISS property. The rele- 
vant definitions are as follows. We say that (10) is integral input-to-state 
stable (iISS) provided there exist 81,82 £ /Coo o,nd (3 £ ICC such that 

rt+t 

8i(\y(t;t ,y ,a)\) < 0(\y o \,t - Q + / 8 2 {\a{r)\)dr (iISS) 

Jt 

everywhere for all measurable essentially bounded functions a : [0, 00) — > U". 
This condition is less restrictive than ISS since e.g. y = — arctan(y) + u is 
iISS but not ISS [3]. An ilSS-LF for (10) (with controls in TJ^) is then 
defined to be a C 1 function V : W n x [0, 00) — > [0, 00) for which there are 
71)72,74 £ fcoo and a positive definite function 73 (i.e., 73 : [0, 00) — > [0, 00) 
is continuous and zero only at zero) such that Conditions 1-2 in Remark 7 
hold everywhere. This is less restrictive than the ISS-LF condition since 73 
need not be of class /Coo. Arguing as in the proof of Theorem 3 up through 
(28) and solving the appropriate constrained minimum problem to get 73 
shows that V = L3 satisfies the iISS Lyapunov function decay condition 
(namely, Condition 2 from Remark 7) for the error dynamics (18) and the 
control set U" using 73(5) = C\(e~ s — l) 2 /2. Therefore, this system is in 
fact iISS, by a slight variant of the proof of the iISS estimate in [3, Theorem 
1]. We leave the details to the reader. 



8. Stability in the Presence of Several Species. Theorem 3 shows 
that the stability of the reference trajectory (6) is robust with respect to 
small perturbations of the dilution rate and initial concentration. To further 
demonstrate the robustness of our results, we next show that the stability of 
(6) is also maintained when the model (4) is augmented to include additional 
species that are being driven to extinction, in the following sense. 

We assume for simplicity that u\ = u 2 = 0. Consider the augmented 
system 

n 

S = D(t)(l-S)-fx(S)x-J2"i(S)yi 
x = x(jjl(S) - D(t)) (33) 
, Vi = yi{vi{S) - D{t)), i = l,...,n, 
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where [i is as in (5) and is continuous and increasing and satisfies i^(0) = 
for i = 1,2, . . . ,n. The variables y% represent the levels of n additional 
species. We choose D and D Q as in (7) and (8), and we assume z^(l) < D Q 
for i = 1,2, ... ,n. (This assumption is, in a sense, natural because one can 
easily check that it ensures that each species concentration converges to 
zero. Indeed, the fact that, for all t > 0, D(t) > D Q > and n(S(t))x(t) + 
127=1 v i{S{t))yi{t) > ensures, in combination with the inequalities ^j(l) < 
D , that there exists an instant T > and a constant c > such that for all 
t > T and for i = 1, 2, . . . , n, ?/j(i) < —cyi(t). This implies that y^'s converge 
to exponentially.) We show that the transformed error 



(z,i,y) := (S + x - S r - x r ,ln(x) -\n(x r ),y) 



(34) 



between any componentwise positive solution (5, x, y) of (33) and the refer- 
ence trajectory (S r , x r , 0, . . . , 0) = (\ — \ cos(t), \ + \ cos(i), 0, . . . , 0) con- 
verges exponentially to the zero vector as t — ► +oo. 

To this end, notice that in the coordinates (34), the system (33) becomes 



-D(t)z(t) 



ma 



i=l 



(a + z - e«)(a + z r (t) - e^W) 



(35) 



Viiyi{S) ~D(t)), i=l,...,n, 



by the same calculations that led to (18). Set Li(£) = exp(£) — £ — 1 where 

= S r = l/2-cos(t)/4 < 1, 



exp(r) := e r . Since e^ r > 1/4, z r = 1, < z 
and < z — 



z + l-e«<2 + z 2 , we deduce that the derivative of 

4m 



L 3 (£,z) :=Li(e) + 
along the trajectories of (35) satisfies 



aD 



-z 2 



L,? < ma 



< ma 



< 



< 



z(e« 



(a + z — e^)(a + z r — e 



4z 2 



16 ^ 



(a + z — e^)(a + z r — e? r ) 
ma(e^ — l) 2 4m 



8m . 
a 

8m 
a 



z 2 



16(a + l)(a + 2 + z 2 ) 



m , 



4Vm \ - 
—=— ^Vi(S)yi 
VaD ^ 
v i=i 



8m 
aD Q 

8m 
aZX, 



^2^i(S)yi 



i=l 



^2vi(S)yi 



i=l 



1 



16m 



16(a + l)(a + 2 + z 



3m 

z 2 H ^> 

a a£>; 



Li=l 



(36) 
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where the second inequality is by (16) with p = \z\, q = |e^ — 1|, and 
d = 4 and the last inequality used the relation J 2 + K 2 > —2JK for real 
values J and K. On the other hand, since z^(l) < D for each i, the form 
of the dynamics for S and the nonnegativity of /i and the u^s along our 
componentwise positive trajectories imply that there exist e > and T > 
such that (i) S(t) < 1 + e for all t > T and (ii) ^(1 + e) < D a for all 
i = 1, 2, . . . , n. We deduce that, for all i = 1,2, ... , n and for all t >T, 

\j t y\ < MS®) - D(t))y? < fa(l + e )-A,)»? < ^ 

where 5 = D Q — max{i/j(l + e) : i = 1,2, ... ,n} > 0. Hence, each tn(t) 
converges exponentially to zero. 

Next notice that along each pair (£(t),z(t)), the function 



16(o+ l)(a + 2 + z 2 ) 
is positive if and only if £ ^ 0. By (36) and (37), the time derivative of 



2/?, where A := — (38) 

ao 

i=i 

along the trajectories of (35) satisfies 



L 4 < -A(£,z) — z 2 + -^ 



o 



^ n{S)yi 



.i=i 



2ASJ2V 2 

i=i 



< 



- A(f , 5) - + ^ £ , 2 (%? - 2A* £ ^ (39) 

° i=l t=l 



< -A(£,5) 2 2 Vy 2 =: -M(£,z,y) 



a a f— f 

i=i 



provided t > T where T is chosen as above. (The second inequality in (39) 
follows because for any nonnegative a^, we get ctfe < (X^iLi a i) 1 ^ 2 which we 
sum and then square to get (XT=i a «) 2 — n2 /Cr=i a i- ^he ^ as ^ inequality 
used Vi(S(t)) < Ui(l + e) < D Q for all t > T and our choice of A. ) 

It is tempting to surmise from (39) and the structure of L4 that L4 is 
a Lyapunov function for (35) since then we could use standard Lyapunov 
function theory to conclude that (£(t), z(t),y(t)) asymptotically converges to 
zero. However, such an argument would not be technically correct, since the 
state space of (35) is not M n+2 (because the original augmented chemostat 
model (33) is only defined for componentwise nonnegative values of the 
state). Instead, we argue as follows (in which we may assume for simplicity 
that the initial time for the augmented error dynamics is zero). 
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For any t > 0, integrating the last inequality of (39) over [0, t] gives 

U (z(t),i(t),y(t)) -U (i(0),e(0),y(0)) 

/"* /- \ (40) 

<-y o M(^(i),m,y(i))di. 

It follows that, for all t > 0, 

£i (£(*)) < L 4 (s(0),C(0),y(0)) . (4i) 

Therefore £(t) = £(t) + £ r (t) is a bounded function. Similarly, z(t) and y(i) 
are bounded. We deduce that £, z, and the components of y are uniformly 
continuous, since their time derivatives (35) are bounded. Reapplying (40) 
therefore implies 

r+oo 

/ M(i(l),z(l),y(l))dl 
Jo 

is finite. It follows from Barbalat's lemma [15, p. 323] and the structure of 
the function M that z(t),y(t)) —> as t — ► +oo. This establishes our 
stability condition for the multi-species model. 

Remark 10. Notice that — L4 is bounded below by a quadratic of the form 
c\(^,z,y)\ 2 along the trajectories of (35), and that L4 is bounded above and 
below by such quadratics along the trajectories as well, since the trajectories 
are bounded. From this fact and (39), one can deduce that the trajectories 
(£(t),z(t),y(i)) converge exponentially to zero. 

9. Simulation. To validate our convergence result, we simulated the dy- 
namics (9) with the initial values x(0) = 2 and S(0) = 1 and the reference 
trajectory x r (t), using the parameters m = 10 and a = \ and t Q = 0. In this 
case, the lower bound on D(t) provided by (8) is D a = 7/3. It follows from 
Remark 9 that the convergence of x(t) to x r (t) is robust to disturbances that 
are valued in [— u, u] 2 for any positive constant u < min{l,D } = 1, in the 
sense of integral input-to-state stability. Moreover, using the estimate (28), 
one easily checks that in this case, estimate (iISS) on p. 14 above holds with 
&2{r) = 2C27"; cf. the proof of [3, Theorem 1]. For our simulation, we took 
the disturbance u\{t) = 0.5e~* on the dilution rate and 1*2(2) = 0. This gave 
the plot of x(t) and x r (t) against time in Figure 3. Our simulation shows 
that the state trajectory x(t) closely tracks the reference trajectory x r (t) 
even in the presence of small disturbances and so validates our findings. 

10. Conclusions. The chemostat model is a useful framework for modeling 
species competing for nutrients. For the case of one species competing for 
one nutrient and a suitable time-varying dilution rate, we proved stability of 
an appropriate reference trajectory. Moreover, we found that the stability 
was maintained even if the model is augmented with other species that are 
being driven to extinction, or if there are disturbances of appropriately small 
magnitude acting on the dilution rate and input nutrient concentration. 
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Figure 3. State Trajectory Component x(t) (dashed) and 
Reference Trajectory x r (t) (solid) for Chemostat 



Appendix. For completeness, we provide the slight variant of the classical 
ISS arguments needed to finish the proof of Theorem 3. Multiplying through 
(32) by e c ' sl and applying the standard "variation of parameters" formula 
to [t Q ,t] B I i ^ V(£(l),z(l)) (by integrating between to > and t > t a ) gives 

V(i(t),z(t)) < e(*°-*) c '«y(f(to),^o)) + C 2 M^^ , (42) 
where we enlarged C2 without relabeling. We deduce that 



z 2 (t) <ln(l + e 



D a -u 



,(to-t)G 



B V(i(t ),z{t )) + C 2 \u\ [toit] 



where L\ is defined in (19). Since e r — 1 — r > \r 2 and ln(l + r) < r for all 
r > 0, we deduce from the formula for V that 



where = e e _1 ~ r+ Do-« r — 1 . 

In particular, VL G /Coo. From (43) and the inequality \J a + 6 < y^a + Vo, 

|£(i)| < \\le^-W*Sl (|(£(to),*(*o))|) + A /2C 2 |«| [to , t] and (44) 



STABILITY IN THE PERTURBED CHEMOSTAT 



19 



The relations z = S-S r + <£ r (e^ - 1) and e a+b - 1 < \ (e 2a - 1) + \ (e 2b - 1) 
give 

\S(t)-S r (t)\ < \~ z (t)\ + (e^-l) 



< yjeto-wtn^n (|(f(to),S(to))|) + ^c 2 ^\ u \ [toA 




The desired ISS estimate (11) now follows immediately from this last 
inequality and (44) with the choices 

p(s, t) = A^n( s )exp(-C 5 t){l + ^} + exp (Vfi(s)exp(-C 5 t)) - 1 
and 7(r) = (1 + (L> D - u)/k) r + exp(4^r) - 1. 

This completes the proof of Theorem 3. 
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